{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.6 DG/HDG splitting\n",
    "\n",
    "When solving unsteady problems with an operator-splitting method it might be benefitial to consider different space discretizations for different operators.\n",
    "\n",
    "For a problem of the form\n",
    "\n",
    "$$\n",
    "  \\partial_t u + A u + C u = 0\n",
    "$$\n",
    "\n",
    "We consider the operator splitting:\n",
    "\n",
    "* 1.Step: Given $u^0$, solve $t^n \\to t^{n+1}$: $\\quad \\partial_t u + C u = 0 \\Rightarrow u^{\\frac12}$\n",
    "* 2.Step: Given $u^{\\frac12}$, solve $t^n \\to t^{n+1}$: $\\quad \\partial_t u + A u = 0 \\Rightarrow u^{1}$\n",
    "\n",
    "This splitting is only first order accurate but allows to take different time discretizations for the substeps 1 and 2. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In this example we consider the Navier-Stokes problem again (cf. [3.2](../unit-3.2-navierstokes/navierstokes.ipynb)) and want to combine \n",
    "\n",
    "* an $H(div)$ -conforming Hybrid DG method (which is a very good discretization for Stokes-type problems) and\n",
    "* a standard **upwind DG** method for the discretization of the convection\n",
    "\n",
    "The weak form is: Find $(\\mathbf{u},p):[0,T] \\to (H_{0,D}^1)^d \\times L^2$, s.t.\n",
    "\\begin{align}\n",
    "\\int_{\\Omega} \\partial_t \\mathbf{u} \\cdot v + \\int_{\\Omega} \\nu \\nabla \\mathbf{u} \\nabla \\mathbf{v} + \\mathbf{u} \\cdot \\nabla \\mathbf{u} \\mathbf{v} - \\int_{\\Omega} \\operatorname{div}(\\mathbf{v}) p &= \\int f \\mathbf{v}  && \\forall \\mathbf{v} \\in (H_{0,D}^1)^d, \\\\ \n",
    "- \\int_{\\Omega} \\operatorname{div}(\\mathbf{u}) q &= 0 && \\forall q \\in L^2, \\\\\n",
    "\\quad \\mathbf{u}(t=0) & = \\mathbf{u}_0\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Again, we consider the benchmark setup from http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html . The geometry:\n",
    "\n",
    "![](http://www.featflow.de/media/dfg_bench2_2d/geometry.png)\n",
    "\n",
    "The viscosity is set to $\\nu = 10^{-3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import netgen.gui\n",
    "%gui tk\n",
    "from math import pi\n",
    "from ngsolve import *\n",
    "from netgen.geom2d import SplineGeometry\n",
    "geo = SplineGeometry()\n",
    "geo.AddRectangle( (0, 0), (2, 0.41), bcs = (\"wall\", \"outlet\", \"wall\", \"inlet\") )\n",
    "geo.AddCircle ( (0.2, 0.2), r=0.05, leftdomain = 0, rightdomain = 1, bc = \"cyl\" )\n",
    "mesh = Mesh( geo.GenerateMesh(maxh = 0.08) )\n",
    "order = 3\n",
    "mesh.Curve(order)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For the HDG formulation we use the product space with\n",
    "\n",
    "* $BDM_k$: $H(div)$ conforming FE space (degree k)\n",
    "* Vector facet space: facet functions of degree k (vector valued and only in tangential direction)\n",
    "* piecewise polynomials up to degree $k-1$ for the pressure\n",
    "\n",
    "![](resources/stokeshdghdiv.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## HDG spaces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "V1 = HDiv ( mesh, order = order, dirichlet = \"wall|cyl|inlet\" )\n",
    "V2 = FESpace ( \"vectorfacet\", mesh, order = order, dirichlet = \"wall|cyl|inlet\" )\n",
    "Q = L2( mesh, order = order-1)\n",
    "V = FESpace ([V1,V2,Q])\n",
    "\n",
    "u, uhat, p = V.TrialFunction()\n",
    "v, vhat, q = V.TestFunction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Stokes discretization / initial conditions\n",
    "The bilinear form to the HDG discretized Stokes problem is:\n",
    "\\begin{align*}\n",
    "  K_h((\\mathbf{u}_T,\\mathbf{u}_F,p),(\\mathbf{v}_T,\\mathbf{v}_F,q) :=\n",
    "    & \\displaystyle \\sum_{T\\in\\mathcal{T}} \\left\\{ \\int_{T} \\nu {\\nabla} {\\mathbf{u}_T} \\! : \\! {\\nabla} {\\mathbf{v}_T} \\ d {x} \\right. \\\\ \n",
    "    & \\qquad \\left. - \\int_{\\partial T} \\nu \\frac{\\partial {\\mathbf{u}_T}}{\\partial {n} }  [ \\mathbf{v} ]_t \\ d {s}\n",
    "    - \\int_{\\partial T} \\nu \\frac{\\partial {\\mathbf{v}_T}}{\\partial {n} }  [ \\mathbf{u} ]_t \\ d {s}\n",
    "      + \\int_{\\partial T} \\nu \\frac{\\lambda k^2}{h} [\\mathbf{u}]_t [\\mathbf{v}]_t \\ d {s}  \\right\\}\\\\\n",
    "          & - \\int_{\\Omega} \\operatorname{div}(\\mathbf{u}) q \\ d {x}\n",
    "           - \\int_{\\Omega} \\operatorname{div}(\\mathbf{v}) p \\ d {x}\n",
    "\\end{align*}\n",
    "where $[\\cdot]_t$ is the tangential projection of the jump $(\\cdot)_T - (\\cdot)_F$.\n",
    "\n",
    "The mass matrix is simply\n",
    "\\begin{align*}\n",
    "  M_h((\\mathbf{u}_T,\\mathbf{u}_F,p),(\\mathbf{v}_T,\\mathbf{v}_F,q) := \\int_{\\Omega} \\mathbf{u}_T \\cdot \\mathbf{v}_T d {x}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "nu = 0.001\n",
    "alpha = 4\n",
    "\n",
    "gradu = CoefficientFunction ( (grad(u),), dims=(2,2) )\n",
    "gradv = CoefficientFunction ( (grad(v),), dims=(2,2) )\n",
    "\n",
    "n = specialcf.normal(mesh.dim)\n",
    "h = specialcf.mesh_size\n",
    "\n",
    "def tang(vec):\n",
    "    return vec - (vec*n)*n\n",
    "\n",
    "a = BilinearForm ( V, symmetric=True)\n",
    "a += SymbolicBFI ( nu*InnerProduct ( gradu, gradv) )\n",
    "a += SymbolicBFI ( nu*InnerProduct ( gradu.trans * n,  tang(vhat-v) ), element_boundary=True )\n",
    "a += SymbolicBFI ( nu*InnerProduct ( gradv.trans * n,  tang(uhat-u) ), element_boundary=True )\n",
    "a += SymbolicBFI ( nu*alpha*order*order/h * InnerProduct ( tang(vhat-v),  tang(uhat-u) ), element_boundary=True )\n",
    "a += SymbolicBFI ( -div(u)*q -div(v)*p )\n",
    "a.Assemble()\n",
    "\n",
    "m = BilinearForm(V , symmetric=True)\n",
    "m += SymbolicBFI( u * v )\n",
    "m.Assemble()\n",
    "\n",
    "f = LinearForm ( V )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As initial condition we solve the Stokes problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(V)\n",
    "\n",
    "U0 = 1.5\n",
    "uin = CoefficientFunction( (U0*4*y*(0.41-y)/(0.41*0.41),0) )\n",
    "gfu.components[0].Set(uin, definedon=mesh.Boundaries(\"inlet\"))\n",
    "\n",
    "\n",
    "invstokes = a.mat.Inverse(V.FreeDofs(), inverse=\"umfpack\")\n",
    "\n",
    "res = f.vec.CreateVector()\n",
    "res.data = f.vec - a.mat*gfu.vec\n",
    "gfu.vec.data += invstokes * res\n",
    "\n",
    "Draw( gfu.components[0], mesh, \"velocity\" )\n",
    "Draw( Norm(gfu.components[0]), mesh, \"absvel(hdiv)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Application of convection\n",
    "\n",
    "In the operator splitted approach we want to apply only operator applications for the convection part.\n",
    "Further, we want to do this in a usual DG setting. \n",
    "As a model problem we use the following procedure:\n",
    "\n",
    "* Given $(\\mathbf{u},p)$ in HDG space: project into $\\hat{\\mathbf{u}}$ in usual DG space\n",
    "* Solve $\\partial_t \\hat{\\mathbf{u}} = C \\hat{\\mathbf{u}}$ by explicit scheme (involves convection evaluations and mass matrix operations only)\n",
    "* Solve Unsteady Stokes step with r.h.s. from convection sub-problem. To this end evaluate mixed mass matrix $\\int \\hat{\\mathbf{u}} \\cdot \\mathbf{v}$ to obtain a functional on the HDG space\n",
    "\n",
    "For the projection steps we use mixed mass matrices:\n",
    "\n",
    "* $M_m$ : $HDG \\times DG \\to \\mathbb{R}$\n",
    "* $M_m^T$ : $DG \\times HDG \\to \\mathbb{R}$\n",
    "* $M_{DG}$ : $DG \\times DG \\to \\mathbb{R}$ (block diagonal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### mixed mass matrices\n",
    "To set up mixed mass matrices we use a bilinear form with two different FESpaces. \n",
    "\n",
    "We do not assemble the matrices as we will only need the matrix-vector applications of $M_m$, $M_m^T$ and $M_{DG}^{-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "VL2 = L2(mesh, dim=mesh.dim, order=order)\n",
    "uL2 = VL2.TrialFunction()\n",
    "vL2 = VL2.TestFunction()    \n",
    "\n",
    "gfuL2 = GridFunction(VL2)\n",
    "\n",
    "bfmixed = BilinearForm ( V, VL2, nonassemble=True )\n",
    "bfmixed += SymbolicBFI ( vL2*u )\n",
    "\n",
    "bfmixedT = BilinearForm ( VL2, V, nonassemble=True)\n",
    "bfmixedT += SymbolicBFI ( uL2*v )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### convection operator\n",
    "* convection operation with standard Upwinding\n",
    "* No set up of the matrix (only interested in operator applications)\n",
    "* for the advection velocity we use the $H(div)$-conforming velocity (which is pointwise divergence free)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "vel = gfu.components[0]\n",
    "convL2 = BilinearForm(VL2, nonassemble=True )\n",
    "convL2 += SymbolicBFI( (-InnerProduct(grad(vL2).trans * vel, uL2.trans)) )\n",
    "un = InnerProduct(vel,n)\n",
    "upwindL2 = IfPos(un, un*uL2, un*uL2.Other(bnd=uin))\n",
    "convL2 += SymbolicBFI( InnerProduct (upwindL2, vL2-vL2.Other()), VOL, skeleton=True )\n",
    "convL2 += SymbolicBFI( InnerProduct (upwindL2, vL2), BND, skeleton=True )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### solution of convection steps\n",
    "We now define the solution of the convection problem for an initial data $u$ in the HDG space. The return value (\"res\") is $M_m \\hat{u}$ where $\\hat{u}$ is the solution of several explicit Euler steps of the convection problem\n",
    "$$\n",
    "  \\partial_t \\hat{u} = - M_{DG}^{-1} C \\hat{u}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def SolveConvectionSteps(gfuvec, res, tau, steps):\n",
    "    bfmixed.Apply (gfuvec, gfuL2.vec) \n",
    "    VL2.SolveM(gfuL2.vec, CoefficientFunction(1))\n",
    "    conv_applied = gfuL2.vec.CreateVector()\n",
    "    for i in range(steps):\n",
    "        convL2.Apply(gfuL2.vec,conv_applied)\n",
    "        VL2.SolveM(conv_applied, CoefficientFunction(1))\n",
    "        gfuL2.vec.data -= tau/steps * conv_applied\n",
    "        #Redraw()    \n",
    "    bfmixedT.Apply (gfuL2.vec, res)\n",
    "    \n",
    "#SolveConvectionSteps(gfu,res,0.01,1)    \n",
    "#Draw(gfuL2, mesh, \"velocity(L2)\")\n",
    "#Draw(Norm(gfuL2), mesh, \"absvel(L2)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Operator splitting method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# initial values again:\n",
    "res.data = f.vec - a.mat*gfu.vec\n",
    "gfu.vec.data += invstokes * res\n",
    "t = 0\n",
    "tend = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  t = 1.000000000000000775"
     ]
    }
   ],
   "source": [
    "tend += 1\n",
    "substeps = 10\n",
    "tau = 0.01\n",
    "\n",
    "mstar = m.mat.CreateMatrix()\n",
    "mstar.AsVector().data = m.mat.AsVector() + tau * a.mat.AsVector()\n",
    "inv = mstar.Inverse(V.FreeDofs(), inverse=\"umfpack\")\n",
    "\n",
    "while t < tend:\n",
    "    SolveConvectionSteps(gfu.vec, res, tau, substeps)\n",
    "    res.data -= mstar * gfu.vec\n",
    "    gfu.vec.data += inv * res\n",
    "    t += tau\n",
    "    print (\"\\r  t =\", t, end=\"\")\n",
    "    Redraw(blocking=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
